#ifndef GPU_PLANE_SPOKE_DART_H
#define GPU_PLANE_SPOKE_DART_H

#include <cuda.h>
#include <curand_kernel.h>
#include "vectorGpu.h"
#include "dataType.h"
#include "gpuDart.h"

class GPUPlaneSpokeDart : public GPUDart
{
public:
	GPUPlaneSpokeDart();
	virtual ~GPUPlaneSpokeDart();

	virtual bool init(unsigned int m, DATA_TYPE radius, unsigned int d);
	virtual void setupDarts();
	virtual void updateIntervals();
	virtual void updateIntervals2();
	virtual void pickPoint();
	virtual void addSamples();
	virtual bool terminationCheck();
	virtual void downloadSamples();
	virtual unsigned int getNumGeneratedSamples();

	unsigned int getNextFront() const { return h_numNextFront_; }
	unsigned int getNumFront() const { return h_numFront_; }

protected:
	GPUMDVector d_directions_;
	GPUMDVector d_computedSamples_;

	unsigned int* d_numSamples_;
	unsigned int* d_numFront_;
	unsigned int* d_failureCount_;
	unsigned int* d_nextFront_;

	unsigned int h_numFront_;
	unsigned int h_numNextFront_;

	// plane spoke loop data
	GPUMDVector d_lastDarts_;
	GPUMDVector d_direction1_;
	GPUMDVector d_direction2_;
};

namespace GPUPlaneSpokeDartCuda
{
	__global__ void initRandomStates(curandState* randStates, unsigned int* pNextFront, GPUDart::GPUMDVector samples, VectorGPU* pFirstSample);
	__global__ void setupGPUDarts(unsigned int* pNumFronts, unsigned int *d_numSamples, GPUDart::GPUMDVector samples, curandState* randStates, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions,
		GPUDart::GPUInterval* intervals, GPUDart::GPUMDVector lastDarts, GPUDart::GPUMDVector directions1, GPUDart::GPUMDVector directions2);
	__global__ void spokeDartLoop(unsigned int* pNumFronts, unsigned int *d_numSamples, GPUDart::GPUMDVector samples, curandState* randStates, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions);
	__global__ void updateIntervalsKernel(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, GPUPlaneSpokeDart::GPUInterval* intervals,
		unsigned int *d_numSamples, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions);
	__global__ void updateIntervals2Kernel(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, GPUPlaneSpokeDart::GPUInterval* intervals,
		unsigned int *d_numSamples, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions);
	__global__ void pickPointKernel(unsigned int* pNumFronts, curandState* randStates, GPUPlaneSpokeDart::GPUInterval* intervals, 
		GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions, GPUDart::GPUMDVector computedSamples);
	__global__ void addSamplesKernel(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, 
		unsigned int *pNumSamples, unsigned int* pFailureCounts, GPUDart::GPUMDVector computedSamples);
	__global__ void updateFront(unsigned int* pNumFronts, unsigned int *pNumSamples, unsigned int* pFailureCounts, unsigned int* pNextFront, unsigned int numBlocks,
		GPUDart::GPUMDVector lastDarts, GPUDart::GPUMDVector directions1, GPUDart::GPUMDVector directions2);

	__inline__ __device__ bool raySphereIntersect(const DATA_TYPE (*sphereCenter)[NUM_THREADS], unsigned int centerIndex,
		DATA_TYPE sphereRadius, const DATA_TYPE* rayOrigin, const DATA_TYPE* rayDirection, 
		DATA_TYPE& t0, DATA_TYPE& t1)
	{
		//Compute A, B and C coefficients
		DATA_TYPE oc[DIMENSIONS];
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			oc[d] = rayOrigin[d] - sphereCenter[d][centerIndex];

		DATA_TYPE a = 0;
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			a += rayDirection[d] * rayDirection[d];

		DATA_TYPE b = 0;
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			b += oc[d] * rayDirection[d];
		b *= 2.0;

		DATA_TYPE c = -sphereRadius * sphereRadius;
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			c += oc[d] * oc[d];

		//Find discriminant
		DATA_TYPE disc = b * b - 4 * a * c;

		// if discriminant is negative there are no real roots, so return 
		// false as ray misses sphere
		if (disc < 0)
			return false;

		// compute q as described above
		DATA_TYPE distSqrt = sqrt(disc);
		DATA_TYPE q;
		if (b < 0)
			q = (-b - distSqrt) / 2;
		else
			q = (-b + distSqrt) / 2;

		// compute t0 and t1
		t0 = q / a;
		t1 = c / q;

		// make sure t0 is smaller than t1
		if (t0 > t1)
		{
			// if t0 is bigger than t1 swap them around
			DATA_TYPE temp = t0;
			t0 = t1;
			t1 = temp;
		}

		return true;
	}
};

#endif